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Abstract 

We simulate the hadroproduction of a tt pair in association with a bb pair at 
14TeV LHC using the PowHel package. We use the generated events, stored 
according to the Les-Houches event format, to make predictions for differential 
distributions formally at the next-to-leading order (NLO) accuracy and we 
compare these to existing predictions accurate at NLO. 
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1 Introduction 



According to the latest, most precise measurements, the t-quark mass is m t = 173.5±0.6± 
0.8 GeV [TJ, indicating that the Yukawa coupling of the t-quark, y t = 0.997 ±0.008 equals 
one with better than 1 % accuracy. This remarkable result suggests that it is important to 
measure this coupling also directly as precisely as possible. After the recent discovery of a 
Higgs-like particle at the LHC (21 13] , the focus of analyses is shifting towards determining 
whether the couplings of this new particle to other fundamental particles agree with the 
Standard Model (SM) expectations. The measurement of the tiH coupling is important 
also because deviation from the SM expectation provides signal of physics beyond the SM. 
If the Higgs-boson is indeed the particle discovered at the LHC recently then its mass is 
approximately 125 GeV. Such a Higgs-boson does not decay into t-quarks, therefore the 
best process to measure the tiH coupling in a model independent way is the study of 
hadroproduction of the Higgs-boson in association with a tt pair [H E]. 

The produced Higgs-boson decays immediately and only its decay products can be 
observed. According to the SM predictions, such a Higgs-boson decays dominantly into 
a bb pair. Thus one can plan to measure the tiH coupling in the pp — > tiH — > ttbb 
process. Unfortunately, this process has a large background from the direct QCD process 
pp — > ttbb. In order to optimize the experimental selection of the tiH events, one needs 
simulation of the QCD background with high precision. 

The most up to date method to simulate LHC processes with high precision is matching 
NLO QCD predictions to shower Monte Carlo programs (SMC). Presently, two formalisms 
are used frequently to match matrix element (ME) calculations at the NLO accuracy to 
parton showers (PS), the MC@NLO [7] and POWHEG [9] methods. Indeed, both 
methods were used to make predictions for tiH hadroproduction [101 [TT], and the two 
predictions were found in agreement in Ref. jJ2j. Here we make one step further and 
describe our implementation of the pp — > ttbb process within the PowHel framework. 

The PowHel framework combines the POWHEG-Box [13J, a flexible computer framework 
implementing the POWHEG method, and the HELAC-NL0 package [14]. The output of 
PowHel is simulated events stored according to the Les Houches accord [15] (LHE). Those 
events can be fed into any shower Monte Carlo (SMC) program for generating events with 
hadrons. Using the same framework we already provided LHE's for several processes at 
the LHC, namely to the production of a tt-pair in association with a hard object, such as 
a jet pjj], or vector bosons [T7] . 

The scope of this letter is to present the simulation of the ttbb LHE's. The process 
under consideration represents an example of unprecedented complexity in PowHel (and 
also among all the processes considered for matched NLO and PS predictions), which raises 
some technical difficulties, worth of detailed description. Phenomenological analyses are 
left for a more detailed publication. 
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Method 



We use the POWHEG formula [9] to estimate the cross section from events with either 
unresolved or resolved first radiation, 



dcr lhe = B($ B )d$ B 

A(<I> B ,^)+d$ rad A($ B ,M$ 



X 



In Eq. (l| B(§b) denotes the NLO-corrected fully differential cross section belonging to 
the underlying Born configuration <3>£ (the integration over the momentum fractions is 
included implicitly), 



(2) 



B{$ B ) = B($ B ) + V($b) + / d$ rad ^($ R ) 
'dx 



+ 



x 



G (B (®b) + G q ($ b ) 



and A(<& B ,P±) is the POWHEG Sudakov form factor that exponentiates the integral of 
the ratio of the real radiation and Born B(<&b) contributions over the radiation 

phase space, 



A($b,Rl) = exp 



B 



(3) 



Using B(<&b), we obtain the cross section at NLO accuracy as an integral over the Born 
phase space, 

o-nlo = / d$ BJ B($ B ) . (4) 



We used the HELAC-NLO package to generate the crossing symmetric matrix elements 
required in the POWHEG-Box as input. In particular, (i) the squared matrix elements for 
the flavour structures of the Born (gg— > tt bb, qq— > tt bb, qq— > tt bb) and real radiation 
emission (qq— > tt bbg, gg— >■ tt bbg, qg— > tt bbq, gq— > tt bbq, qq— > tt bbg, qg— > tt bbq, 
gq— > tt bbq) subprocesses (qe{u,d,c,s} - we neglect the contribution of b-quarks in the 
initial state), (ii) the colour-correlated and spin-correlated squared matrix elements for 
the Born flavour structures, and (iii) the finite part of the virtual correction contributions 
in dimensional regularization, on the basis of the OPP method [18] complemented by 
Feynman rules for the computation of the QCD R2 rational terms [19]. With this input 
POWHEG-Box generates events with either unresolved or resolved first radiation. Then, one 
can use any shower Monte Carlo (SMC) program for generating events with hadrons. We 
leave the latter together with phenomenological analyses to a separate publication. 
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2.1 



Checks 



The consistency among real-emission, Born, colour-correlated and spin-correlated matrix 
elements was checked in randomly chosen phase space regions by taking the soft- and 
collinear limits of the real-emission squared matrix elements in all possible kinematically 
degenerate channels. 

The process presented here was studied at the NLO accuracy in the literature [20] l2Tj . 
which enabled us to make detailed checks of our calculation. We reproduced all figures of 
Ref. [20] and found complete agreement. 

2.2 Generation cuts and suppression 

As we treat the b-quarks massless, the fully differential Born cross section becomes sin- 
gular if any of the transverse momenta of the b-quarks p_i_,b, or the invariant mass of 
the bb-pair m h ^ vanishes. To make the integral of B finite over the whole Born phase 
space, we introduce generation cuts by requiring pj_ h > p^f' '* and m h ^ > , where 
p(g-c) _ m fe_ c -) — 2 GeV. We checked that these generation cuts are sufficiently low, so 
that the physical predictions are independent of those when selection cuts characteristic 
to physical analyses are superimposed. However, with such low generation cuts the event 
generation is very inefficient because most of the events are generated in regions of the 
phase space where B is large, but those regions are usually not selected in the analyses. 
Therefore, we also introduce a suppression factor [22] for B of the form 



where pf pp) = m^ pp) = 30 GeV. 
2.3 Size of the remnant 

Depending on the particular process being considered in the POWHEG method, the size 
of the R/B ratio in the Sudakov exponent can become much larger than its collinear 
or soft approximation. In such cases the generation of the radiation becomes highly 
inefficient. This is especially true for the process under consideration. In order to speed- 
up calculations the POWHEG-Box contains the option of separating the real emission phase 
space into two disjoint regions, one containing the singular regions while in the other the 
squared matrix element remains finite. The latter is called the remnant contribution. 
Thus with this separation the real emission part can be cast into a sum R = R s + R T , 
where R s and R r are the singular and finite contributions, respectively. In the original 
POWHEG-Box implementation this separation is controlled by a step function in the form 



supp 



( 




(5) 
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Rs — (1 — T)R and R T = J-R, where 



j- = \ I, R > X max(C co n , C soft ) . . 

1 , otherwise 

with C co u and C so ft being the sums of the collinear and soft counterterms and x is a f ree 
parameter that can be tuned to control the size of the phase space of the remnant. The 
larger Xi the smaller this region, and disappears completely in the limit x ~~>* 00 • With 
this separation the POWHEG formula in Eq. (Jll) becomes a sum of two terms, 

CLHE = CLHE,s + CLHE,r • (7) 

The first one, <7lhe,s is defined as in Eqs. ([T|-(|3]) with the change that R is replaced with 
Rs everywhere. The second one contains the contribution of the remnant, 



a L HE,r = y d$ Rj R r ($ R ) . (8) 

Similarly, the cross section at NLO accuracy in Eq. Q becomes a sum of two terms, 

C"NLO = C"NLO,s + 0"NLO,r , (9) 

where the first term is obtained as in Eq. Q with the R(&r) — > i? s ($j?) substitution and 
the second term is equal to o"LHE,r- 

In principle the size of the remnant can be controlled freely by the parameter x, but 
practical considerations limit this freedom. Unfortunately, when using the POWHEG 
method, the efficiency of event generation decreases with increasing x because over a 
large part of the phase space the value of the fully differential cross section is much 
smaller than its peak value. This decrease in efficiency depends on the process under 
consideration. For the complex final state of tt bb-production, the size of the phase space 
of the remnant is large if for instance, the value of x is chosen to its default value in the 
POWHEG-Box, x = 5- 

The separation of the real emission contribution does not influence the NLO cross 
section. We show the result of an example computation in Table [TJ These computations 
were done for proton-proton collisions at yfs = 8 TeV, with CT10 parton distribution 
functions (PDF), equal renormalization and factorization scales jd,® — ^Ur, = /^p = mt — 
173.2 GeV, and the virtual part was replaced with the half of the Born squard matrix 
element V — > B/2 (which we call 'fake virtual'). We imposed the following set of cuts: 



1. A track was considered as a possible jet constituent if |r/| < 5, t-quarks were 
excluded from the set of possible tracks. 

2. Jets were reconstructed with the /cj_-algorithm using p™™ = 20 GeV and R = 0.4. 

3. We required at least two jets, one b- and one 6-jet, with 1 2/t. (5) I < 2.5. 

4. Events with invariant mass of the 66-jet pair below m^ in = 20 GeV were discarded. 
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(Tmt n ffbl 

^ IN J_jVJ L a 
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5 


446 ± 1 


1109 ± 9 


1555 ± 9 


22% 


10 


514 ±2 


1040 ± 9 


1554 ± 10 


17% 


50 


693 ±2 


855 ±7 


1549 ± 7 


1.9% 


oo 


1577 ±16 


0±0 


1577 ± 16 


~0% 



Table 1: The two contributions to ctnlo according to the decomposition of Eq. ([9]) 
computed with fake virtual. The last column shows the ratio of LHE's without resolved 
radiation to those with resolved first radiation. 



The separation of the real emission into the singular and remnant contributions however, 
does influence the distributions computed from the generated LHE's. We study this 
dependence in detail in the next section. Here we only point out that the ratio of the 
events with Born kinematics to those with first radiation depends on the value of %. In 
the last column of Table [l] we show this ratio, and find that for the default separation in 
the POWHEG-Box, \ = 5, the number of generated events with Born kinematics is too large, 
making the transverse momentum distribution of the hardest non-6 jet (rljJ incorrect. 
Increasing x, this ratio decreases and becomes approximately % in the limit x ~~ °°i 
and the p±j 1 -distribution fits close to the NLO prediction for large transverse momenta. 
(At small p±j 1 we see the effect of the Sudakov-damping.) This suggests that one should 
use at least x = 50 to separate the remnant without introducing too large fraction of 
events with Born kinematics and, as a result, incorrect distribution of the hardest non-6 
jet. 

2.4 Precision and efficient evaluation of loop amplitudes 

In the computations above we used a fake virtual part V = B/2. The reason for such a fake 
computation is the high degree of difficulty of evaluating the loop amplitudes. The high- 
rank of the tensor loop-integrals makes the numerical computation of the loop amplitude 
unstable, and thus, unreliable when double precision arithmetics is used in CutTools [T5] , 
as implemented in HELAC-lloop (part of HELAC-NL0). In order to control numerical insta- 
bilities we employ an M = M test as implemented in CutTools. For a given numerator 
we determine the scalar-integral coefficients using double precision arithmetics. We check 
the accuracy of the integrand by reconstructing it using all the coefficients (and spurious 
terms) with a randomly chosen loop momentum. If the reached relative accuracy is worse 
than 10~ 4 , we pass the same phase space point in double-double precision (computed si- 
multaneously with the double-precision version) to HELAC-lloop@dd to recalculate all the 
coefficients (using the multiple precision part of CutTools). The HELAC-lloop@dd code is 
a straightforward extension of HELAC-lloop to double-double precision using the package 
QD [23]. This way we solved all numerical instabilities in the computation of the virtual 
correction. The above procedure however, makes the numerical computation of the loop 
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Figure 1: Transverse momentum distribution of the leading non-6 jet at the LHC at 8 TeV 
with fake virtual contribution at three different value of the parameter x that controls 
the size of the remnant contribution (dashed: x — 5, long dashed: x = 10, dot-dashed: 
X = 50), compared to the NLO prediction (solid) with the same fake virtual part. The 
inset shows the same plot with range restricted to small transverse momenta. 



amplitudes rather cumbersome, requiring much CPU time. This problem is magnified by 
the hit-and-miss procedure, the way event generation is performed in the POWHEG-Box. 
During NLO integration the inclusive NLO cross section is computed simultaneously with 
the maximal value of B 

S max = maxS($ B ). (10) 

When the underlying Born kinematics is generated for an event, it is accepted if 
B(&b) > £,B mSLX , where £ is a random number picked between zero and one. Using 
this hit-and-miss method, the B function has to be evaluated several times to find a 
suitable phase space point which is selected. If B is calculated multiple times per event, 
event generation becomes highly inefficient because B contains the virtual part, whose 
evaluation is very time consuming. To improve the efficiency, we first generate events 
with the fake virtual. We select those events that pass the selection cuts, whose weights 
in a final step are reweighted by the true virtual contribution. This way we can reduce 
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Figure 2: Distribution of weights before and after reweighting. 



the computation of the virtual contributions to few hundred thousand phase space points. 
The plots shown in the next section contain 200 k points. 

The POWHEG method produces unweighted events. When we generate events with 
the fake virtual part, unweighting must be done using the NLO cross section obtained 
also with the fake virtual. After the events that pass the selection cuts are collected, we 
reweight those 



rv \ V=B/2 vv 



w 



V=B/2 



B \v=b/^b) 



11) 



where VV| V/=B ^ 2 is the weight with fake virtual, and W is the weight with the true virtual 
part. Clearly, the new weight becomes dependent upon the underlying Born kinematics. 

The majority of events generated by the POWHEG-Box have either positive or negative, 
but in magnitude equal weights. As seen on Fig. [2} the relat ive fraction of the negative 
weights is at the percent level. (Due to the suppression in Eq. (|5j there are also several 
large weight events starting at about 50 in Fig.[2j) The drawback of the reweighting 
procedure is the loss of equal weights: although about 90 % of the events have equal 
positive weights, LHE's with many different weights also appear. 
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3 Predictions from the LHE's 



Using PowHel one can make predictions at five different stages in the evolution of the 
final state: (i) at the parton level using NLO accuracy, (ii) from the LHE's, formally 
at the NLO accuracy, (iii) after decay of the heavy particles, (iv) after parton shower, 
(v) at the hadron level after full SMC. While in an experimental analysis the last option 
is the most useful, for examining the effect of the different stages in the evolution, and 
also for checking purposes, it is useful to study the predictions at intermediate stages. In 
the previous section we mentioned that we made comparision with existing predictions 
at stage (i) and found agreement. We now show predictions at stage (ii), and discuss 
how these depend on the separation of the remnant (on x) an d on the choice of the scale 
fio, and leave the phenomeno logical analyses at the remaining stages to a more detailed 
publication. 

Our selection cuts are those used in setup 1 of Ref. [21], which amounts to the cuts 
(1-3) above supplemented with an invariant mass cut of the bb-jet pair, m h i > 100 GeV. 
In Figs.|3]-[5ja we show our predictions from the LHE's as bands obtained by varying the 
scale in the range [/zo/2,2/x ] around the fixed default scale fiQ = m t and the predictions 
at NLO accuracy with scale varied in the same range. The LHE's were obtained with 
remnant separation parameter \ — 5, but for comparison we also show the predictions 
with the default scale from LHE's generated with x — 50, which describes the distribution 
of the leading non-b jet more accurately. In Figs.[3]-[5jb we show our predictions from the 
LHE's as bands obtained by varying the scale in the range [/xo/2, 2/xo] around a dynamical 
scale (Mq = m t ^p± ! bP±j ) - The LHE's were obtained with x — 50. For comparison on each 
plot we also show the NLO predictions of Ref. [21] obtained with their default dynamical 
scale, and the band of our NLO predictions obtained by varying the fixed scale in the 
range [/i /2, 2/x ]. 

We present kinematical distributions of the 6-jet pair: transverse momentum in Fig.|3j 
rapidity in Fig.|4j and invariant mass in Fig.[5j Our first observation is that the predictions 
from the LHE's show much smaller scale dependence than the corresponding predictions 
at NLO accuracy. Nevertheless, the latter bands overlap with the former ones. 

Looking at the distributions obtained with the fixed scale, it is remarkable that the 
predictions from the LHE's obtained with x — 5 coincide almost exactly with the NLO 
predictions of Ref. [21] obtained with a dynamical scale, which is true for all kinematical 
distributions presented in Ref. [21]. This suggests that these LHE's can be used if one is 
interested in distributions related to observables of particles emerging from the primary b- 
quarks. Unfortunately, as seen in Fig.[T]the events obtained with x — 5 do not describe the 
non-b radiation. Instead, p± ! j 1 is described correctly if x — 50. However, the predictions 
from the LHE's obtained with x — 50 and fixed scale fall below the distributions obtained 
with x — 5, though in agreement with it within scale uncertainty. There is one exception, 
the invariant mass distribution of the 66-jet pair in the range of [100,150] GeV, relevant 
to Higgs studies. It is encouraging that in this range the various predictions are robust 
against variations of the parameters. 
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Figure 3: Transverse momentum distribution of the 66-jet pair at the LHC at y/s = 14 TeV 
using PowHel. Distributions from LHE's are denoted PowHel LHE, while those at NLO 
accuracy by PowHel NLO using a fixed scale hq = m t . The shaded band corresponds to 
cross section obtained with varying the scale in the range [/to/2, 2/xq]- F° r comparison we 
also show the NLO predictions from Ref. (BDDP NLO). (a) PowHel LHE's generated 
using a fixed scale fiQ = m t and with remnant separation x = 5 (solid line with solid 
band as scale dependence), and with x = 50 (dash-dotted), (b) PowHel LHE's generated 
using a dynamical scale //jj = rn t ^p± t bP±,bj an d with x = 50 (solid line with solid band as 
scale dependence). The left scale on the lower panels show the ratio of the bands of scale 
dependences to the distribution of the PowHel LHE's while the right scale on the lower 
panels shows the ratio of the NLO prediction to that obtained from the same PowHel 
LHE's generated with (a) fixed scale and x — 5, (b) dynamical scale and x — 50. 
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Figure 4: Same as Fig.[3]for the rapidity distribution of the bb-jet pair. 



It is also convincing that looking at the distributions obtained from the LHE's gener- 
ated with the default dynamical scale and remnant separation parameter x — 50, shown in 
Figs.|3]-[5]b, the predictions are close (within 10%) to the NLO ones of Ref. [2T], obtained 
also with dynamical scale. This observation is valid for the comparison of all kinematical 
distributions from our LHE's and the NLO predictions of Ref. [21] . Thus we conclude 
that these events can be used to produce distributions that agree with the NLO predic- 
tions up to corrections beyond NLO accuracy, and can be used as input in the SMC's for 
generating events at stages (iii-v). 
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Figure 5: Same as Fig.|3|for the invariant mass distribution of the 66-jet pair. 

4 Conclusions 

We have implemented the hadroproduction of a tt-pair in association with a bb-pair, 
pp — > ttbb, in the PowHel framework, which can be used to generate LHE files. We 
discussed the technical subtleties of such event simulations, in particular, issues related 
to the efficient generation of the events for such a complex process: (i) the advantage of 
separating the phase space into singular and remnant regions and the dependence of the 
predictions on this separation, (ii) the implementation of higher precision arithmetics in 
computing the loop amplitudes with sufficient precision, (iii) the advantage of performing 
the event generation with fake virtual part and the method of reweighting. 

The event files produced by PowHel, together with further detail and results of our 
project, the corresponding version of the program are available at |http: //grid. kfki.hu/twiki 
|/bin / view / DbTheory | We are confident that these LHE's are suitable input into SMC's 
to produce distributions at the hadron level needed for experimental analyses. Thus us- 
ing those events one can make a detailed analysis of tiH signal and ttbb background 
predictions at the hadron level. 
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